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ABSTRACT 


This  paper  discusses  a number  of  general  concepts  relating  to 
relaxation  procedures,  methods  of  “steepest  descent,”  an<l  itera- 
tion  procedures.  The  relationship  between  these  is  indicated 
for  certain  systematic  ways  in  which  relaxation  patterns  may  be 
developed.  An  extrapolation  procedure  is  developed  for  prob- 
lems in  which  the  successive  approximations  technique,  in  its 
usual  form,  diverges.  The  procedure  is  a generalization  of  an 
observation  by  Hartree  and  others  that  in  many  cases  such  di- 
vergent problems  can  be  treated  by  considering  the  divergence 
to  be  a geometric  one.  The  study  of  upper  and  lower  bounds  to 
the  errors  in  relaxation  procedures  is  a second  major  part  of  the 
paper.  General  procedures  arc  developed  for  determining  the 
maximum  and  minimum  errors  in  any  of  the  quantities  for  which 
bounds  to  an  influence  function  can  be  determined.  Use  of  rhe 
theorems  derived  makes  it  possible,  for  example,  to  estimate  the 
error  due  to  all  sources  including  “round  off”  at  any  stage  in  the 
numerical  solution  of  Laplace’s  or  Poisson’s  equation 


INTRODUCTION 


In  this  equation,  x.  may  be  considered  as  the  deflection 
at  a point  i,  —6,  as  the  load  at  the  point,  and  r . as  a resid- 
ual which  becomes  zero  if  all  the  values  of  x are  chosen 
correctly. 

It  is  convenient  to  deai  with  a correction  pattern  of  de- 
flections Uj  which  produces  a change  in  the  residuals  of 
magnitude  p,-, 

Yjaikuk  = pi  t2] 

The  quantity  x + u may  be  a better  approximation,  in 
general,  to  the  solution  than  is  x.  A better  approximation 
can  always  be  obtained  with  the  quantity 

x + tu  [3] 

which  involves  the  parameter  or  multiplier  t.  The  solution 
of  the  equations  corresponding  to  r = 0 is  given  by  the 
addition  of  h to  x,  and  is  defined  as  <f>  in  the  relation 


1.  Summary 

This  paper  describes  certain  formal  procedures  for  the 
numerical  solution  of  problems  which  can  be  stated  in  the 
form  of  a set  of  simultaneous  linear  algebraic  equations, 
although  some  of  the  applications  extend  beyond  these 
limitations.  Consideration  is  given  to  a method  for  inves- 
tigating upper  and  lower  bounds  to  the  numerical  values 
of  the  solution  obtained  at  any  stage.  The  closeness  of 
the  bounds  depends  upon  the  choice  of  a correction  pat- 
tern, but  general  rules  for  selecting  the  best  pattern  can- 
not yet  be  given. 

An  investigation  is  made  of  the  close  relationship  be- 
tween extrapolation  processes  which  can  be  used  to  speed 
up  the  convergence  of  iterative  procedures  and  "methods 
of  steepest  descent"  which  are  used  in  formal  relaxation 
procedures.  Although  the  convergence  of  these  procedures 
is  slower  than  for  methods  which  involve  the  continual  ex- 
ercise of  judgment,  it  seems  possible  to  code  the  formal 
methods  for  use  with  automatic  digital  computers. 

2.  General  Concepts,  Definitions,  and  Notation 

It  is  assumed  that  the  reader  is  familiar  with  the  more 
common  details  of  relaxation  and  iteration  procedures. 

'X'e  consider  the  set  of  equations  of  which  the  following  is 

a particular  equation: 


(f> . = x + h , [4] 

1 1 1 

The  ordinary  iteration  procedure  is  essentially  the 
determination  of  a new  set  of  values  x from  a given  set 
x by  means  of  the  following  scheme.  In  the  i th  equation 
of  Ll],  set  rt-  = 0,  and  take  all  x ^ as  given  except  for 
k = i.  Choose  x ■ so  as  to  satisfy  the  equation.  Repeat 
for  each  equation  in  [ 1 j.  In  general,  x i will  be  different 
from  x,-;  in  fact 

x'i  = xi  ~ r/aii 

This  can  be  rewritten  as  follows: 


x 


in  which 


"  1 2 * *  5ik>k  + bi 


8 . , = 0 for  i 4s  k, 

ik 

8ik  = 1 for  i =-  k. 


[51 
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We  shall  choose  to  write  the  relation 


X ' *=  X(  r Ut  (6] 

to  indicate  the  change  in  x produced  by  the  single  step  in 
the  iteration  process.  Similarly,  if  we  start  with  a value 
y instead  of  x,  we  obtain  a new  value 

y ' - y * v . [7] 

Other  notation  will  be  defined  when  i:  is  introduced. 

UPPER  AND  LOWER  BOUNDS 

3.  Ccnctul  Relations 

The  solution  of  ( 1 1,  with  r « 0,  can  be  stated  as 


<5  **  x.  + h 


(8] 


Now,  suppose  some  quantity  O is  desired  which  is  a 
linear  function  of  <p,  as  follows: 


O -T 

- n 


( / ni  ’i 


(91 


Then,  let  us  define  a new  quantity  p,  such  that 
1 * k 


O 


2 1101 


in  which 


^ i ani7'  1* 


11) 


The  quantities  can  be  considered  as  an  influence 
function  for  the  quantity  0 for  unit  values  of  the  loaa 

In  the  same  way  that  9 corresponds  to  the  true  solution 
d>,  let  us  define  a quantity  .M  corresponding  to  the  approxi- 
mate solution  x, 


\—5X(h-v 


This  follows  from  the  fact  that  the  deflections  x corre- 
sponds to  loads  — A,  ■»  r(  in  the  same  way  that  the  deflec- 
tions 6,  correspond  10  loads  — A,. 

It  follows  immediately  that  the  difference  between  O 
and  '!  can  ! e expressed  as 

Sh-r>n  - XX*'*  (13) 


If  this  difference  is  positive,  then  0 is  definitely  less 
than  M,  and  if  it  is  negative,  O is  greater  than  Al.  It 
follows,  therefore,  that  if  we  know  the  values  of  /),  then 
for  any  set  of  values  x which  implies  a set  of  values  of 
r,  we  could  determine  whether  M is  greater  than  or  less 
than  O.  But  we  do  not  know  fi.  However,  under  certain 
conditions  we  can  determine  bounds  or  upper  ard  lower 
limits  to  fi  readily  enough. 

Let  us  designate  those  values  which  are  upper  bounds 

by  the  symbol  fia  and  lower  bounds  by  (i^ . That  is, 
whether  fi  is  positive  or  negative,  the  relation  holds,  that 

{ 14] 


Next,  we  separate  from  r ^ those  values  which  are  positive 
from  those  which  are  negative  and  designate  thaia  respec- 
tively by  and  r”  . In  investigating  the  sign  of  the  right 

hand  side  of  tl3l>  we  must  use  separately  the  upper  and 
lower  limits  of  B,  as  follows- 

If 


k 

Efc  'TCi 

then 

At  > 0 and  M - At  " 

n — ' n n n 

and  if 

i 

vbnk  \)<° 

then 

At  <0  and  M = .M  b. 

n ~ '«  n » 


(131 


(161 


If  these  inequalities  are  not  satisfied,  the  criterion  fails. 
4.  Use  of  Relaxation  Patterns 

In  the  study  of  bounds  it  will  rarely  occur  that  the 
criteria  of  equations  ( 15)  and  ( 16)  will  be  satisfied  for  a 
paiticular  set  of  residuals  and  given  values  of  upper  and 
lower  limits  to  fi.  We  shall,  therefore,  consider  ways  of 
adding  a proportion  of  a correction  pattern  u to  x in  order 
to  arrive  at  a useful  result.  The  pattern  u corresponds  to 
changes  in  residuals  />,  in  accordance  with  equation  (2l. 
We  seek  a multiplier  ! such  that  the  solution  x + tu  corre- 
sponds to  residuals  t t/f  which  give  us  an  upper  or  lower 
bound  for  O,  according  to(l5l  nrll6). 

For  any  range  in  the  values  of  t for  which  r ^ + tp^  does 

not  change  in  sign,  we  can  derive  from  equation  ( 1 5 1 rhe 
relation 

* k 

X ,rk + v*)’ ’ X^C  - 'pp~ (l?i 


The  quantities  in  parentheses  must  be  ordered  according 
to  the  sign  of  their  sum,  as  indicated.  Now,  those  ”alues 

of  r and  p which  come  from  i >p)*  arc  designated  r ■* 


10 


and  p ’ - even  ti-^ugh  some  of  these  values  may  be  nega- 
tive; and  similarly  for  r^and  ft ' which  come  from 
(r  + tp)  Then  we  can  write: 


E[a\  ».'■'] 

- ±[&  ■v*’**:  '.i  i 


r + tp  .0 

i r t — 

Now,  let  us  define  / as  follows 


/ r=  — r 

; i 


[221 


(231 


from  which  we  determine  a minimum  value  of  /: 


[13] 


Efc 


Ue  must  consider  separately  positive  and  negative  values 
of  p^.  Ve  shall  designate'  the  greatest  algebraic  value  of 

t.  for  positive  values  of  p as  and  for  negative  values 

of  p as  Ta_  . Similarly  we  shall  designate  the  minimum 

algebraic  values  of  for  positive  />  as  7~^  and  for  nega- 
rfc 

tive  p as  l 

Then,  in  order  that  equation  [22]  be  satisfied  for  all 
value::  of  r,  the  following  condition  must  apply; 

rSor"  [24] 


(l9i 


If  the  value  of  I in  [l9l  requires  changes  in  the  sign  of 
any.  + tp^,  the  correspondirg  cnange  in  the  summations 

must  be  made.  However,  by  a ivncess  involving  several 
trials  a value  of  t may  be  obtained  which  corresponis 
therefore  to  the  desired  upper  bound  Al^*3. 

A similar  expression  can  readil)  be  written  for  the 
value  of  i,  from  equation  [ 1 6] , corresponding  to  a lower 

bound  . It  is  clear  that  the  values  obtained  depend 

on  the  patterns  used,  and  in  many  cases  it  may  not  be  j>os- 
sible  to  satisfy  the  criteria  for  a given  pattern. 

In  certain  cases  it  is  convenient  to  use  as  a correction 
pattern  a set  of  values  proportional  to  r.  This  does  net 
require  special  consideration  beyond  the  remark  that  when 
u = x,  the  changes  in  residuals  p arc  given  simply  by  the 
equation 

-b,.  [20] 


and  it  is  only  under  these  restricted  conditions  that  this 
relation  can  be  satisfied. 

Similarly,  in  order  that  the  following  relation  hold: 

r * tp  .<  0 1 2?>1 

i ■ i — 

the  required  condition  is, 


7 \ > t > t 


[26] 


Another  special  case  of  considerable  interest  con- 
cerns bounds  for  the  volume  under  a surface  the  ordinates 
to  which  must  satisfy  given  conditions.  Numerically,  the 
volume  can  be  expressed  as  a set  of  multipliers,  essen- 
tially nearly  equal,  times  the  individual  ordinates,  or  as 
a consequence,  as  a constant  times  the  sum  of  ihe  oiJi- 
nates.  In  other  words,  the  values  of  a in  equations  [9l 
and  [ 1 1 ] are  unity. 

It  follows,  therefore  that 


- s «,■ 


£*, 


[27] 


V 


But  in  cases  where  x nr  u is  a deflection,  and  r or  p is  a 
ioad,  we  can  write  the  relation: 


5.  Special  Cases 

Certain  special  cases  are  of  particular  interest  because 
of  their  simplicity  and  general  application.  For  example, 
in  such  problems  as  the  flexure  of  certain  plates  or  the  de- 
flection of  membranes  it  is  obvious  that  a ioad  atanypoinr 
produces  a deflection  a».  every  point  in  the  same  general 
direction.  For  deflections,  at  least.it  is  clear  that  for  any 
point  n,  a lower  limit  to  fi  is  a positive  number  approach- 
ing zero,  and  an  upper  limit  is  some  positive  number  B.  It 
follows  directly  from  equations  [15]  and  [16]  that: 

‘.I  - A!"  if  all  residuals  arc  positive, 

and  [21] 

A!  »»  A!5 &  if  alt  residuals  are  negative. 

These  conditions  are  not  nearly  as  complicated  as  those 
that  might  be  used  toobtain  more  accurate  bounds,  but  they 
are  extremely  useful  as  a guide  tojudgment  inrelaxation  so- 
lutions for  problems  where  the  erferia  ore  applicable. 

Allien  all  the  influence  values  u.e  positive,  a simpler 
condition  than  [l9i  can  be  derived.  We  seek  a multiplier 
i so  as  to  make,  for  example,  every  value  of  rota!  residual 
{KJSitive,  i.e. 


A V 

\ = - £ r* 


[28] 


- 1 


(.where  u * means  the  deflection  at  i due  to  a unit  load 

at  k).  Te  can  regard  the  differences  between  the  true  de- 
flect i<'o°  and  x as  caused  by  the  residuals  r. 

However,  because  of  Maxwell's  theorem  of  reciprocal 
deflections,  [?8j  can  be  rewritten  as  follows: 


k r = 1 


Then  we  can  derive  the  result: 


£*. 


-££ r*  “* 


- 1 


£ r*  £ “* 


[29l 


[30] 


But  the  last  summation  is  merely  the  deflection  at  k due 
to  unit  loads  at  all  points  i,  and  wc  can  designate  rhis 
by  the  symbol  u.. 

* i 

As  a result  we  have,  if  we  designate  x.  b>  the 
symbol 


11 


* - <?  - E 

where  u,  is  the  deflection  jxitt 

load  of  unity  at  every  point.  But  this  is  precisely  the 
same  equation  as  [lj]  with  /3^  replaced  by  u conse- 
quently the  same  comments  apply,  and  the  criteria  ex- 
pressed by  equations  1 1 5 1 and  [Hi]  can  be  used  to  de- 
termine whether  M is  an  upper  or  lower  bound,  if  we  can 
obtain  bounds  for  u.  In  general  this  is  net  particularly 
difficult  to  do  for  such  problems  as  are  encountered  with 
membranes  and  plates. 

Since  the  determination  of  the  torsion  constant  for  a 
prismatic  bar  requires  the  calculation  of  a volume,  bounds 
- to  the  torsion  constant  can  usually  be  determined  quite 
accurately  from  equations  (3 1 1.  ( 1 5l  and  [16].  Thesegive 
a much  closer  limit  than  can  be  obtained  from  the  condi- 
tion requiring  all  residuals  to  be  of  the  same  sign. 

6.  Hounds  for  Influence  Function 

In  the  use  of  the  procedure  described  herein  it  is  nec- 
essary to  estimate  upper  and  lower  bounds  for  the  influ- 
ence functions  in  the  equations.  The  simplest  influence 
function  is  that  for  the  deflection  at  a particular  point  t:. 
for  which  the  values  of  a in  equations  [$)]  and  [ll]  are 

zero  except  tor  i = n.  This  influence  function  can  be 
obtained  qualitatively  in  most  physical  problems  quite 
readily  as  the  deflection  function  for  a unit  load  at  point 
n.  However,  no  simple  rule  for  the  determination  of  the 
influence  function  can  be  given  for  problems  involving 
arbitrary  systems  of  equations. 

In  physical  problems  which  deal  with  loads  and  de- 
flections, use  can  be  made  of  a simple  theorem  tc  restate 
the  problem  of  determining  an  influence  function  in  terms 
ci  the  problem  of  finding  deflections  for  given  loads. 

This  theorem  is  derived  in  a previous  paper  by  the  author1 
an  1 is  stated  as  follows: 

Where  a particular  effect  O is  a linear  function  of  the 

* ' n 

I 

d»f  lections  r.  at  points  i,  as  in  theequationO  - > A x , 
* ■ r.  i—t  i,  i i 


Uk  '* 


(31] 

ern  produced  bv  a uniform 


bounds  for  tne  influence  function  either  b)  use  of  judg- 
ment or  by  a rough  calculation  and  several  cycles  of  re- 
laxation. It  appears  offhand  to  be  necessary  to  go 
through  the  problem  several  times,  but  this  is  hardly 
warranted,  since  the  trouble  involved  in  closely  delimit- 
ing the  bounds  is  great  enough  to  warrant  instead  the 
calculation  of  more  accurate  figures  with  additional  cy- 
cles of  relaxation  or  iteration. 

EXTRAPOLATION  WITH  ITERATIVE 
PROCEDURES 

7.  Genera!  Relations 

A brief  but  interesting  and  valuable  discussion  of 
•iterative  procedures  has  been  given  by  llartree*  in  which 
he  states  a process  of  improving  or  of  essentially  has- 
tening the  convergence  of  a sequence  of  iterative  cal- 
culations of  numbers.  There  is  described  here  a gen- 
eralization of  this  procedure  which  applies  to  sets  of 
numbers,  and  particularly  to  iterative  solutions  of  simul- 
taneous equations,  although  it  is  applicable  to  other 
pro!  lems. 

Consider  a set  of  numbers  x w-hich  after  one  itera- 

t 

five  cvcle  leads  to  a new-  set  x *?  x + u . A similar 

• i i i 

set  y leads  to  y s-  y + r 

i t - i i . 

We  seek  a set  of  numbers  z leading  to  z'  which  is  to 

i n i 

be  closer  in  general  ro  than  x'  or  y'  is  to  x or  y. 

tthen  the  system  is  linear  or  even  approximately  linear, 

we  proceed  as  ft’ lows: 

z - x + t (y  - x ) 

i i ; 

T - X 4 t (y  - X ) rt  z 

t l It  l 

Now  let  us  choose  to  minimize 

£ [“, + 1 <i , 

This  leads  to  the  condition 


tu  W(i  — u ) [ s2 1 

i i i 

the  expression 

-O]’  1 33 1 


and  if  deflections  and  loads  arc  linearly  related,  the  in- 
fluence on  0 of  a unit  load  a:  any  point  k is  obtained  as 

the  deflection  at  the  point  k due  to  loads  of  magnitude 

A applied  at  the  points  t . 
m 1 1 


l r-  - 


where  u 

t 


i - u 
t i 


[341 


In  problems  to  which  this  theorem  is  applicable,  the 
bounds  for  the  influence  functions  can  be  obtained  in  the 
same  way  as  bounds  for  any  other  problem. 

Consider,  for  example,  the  problem  of  determining 
bounds  to  the  moment  in  a plate  at  a specific  point  due 
to  a given  distribution  of  loading  on  the  plate.  In  solving 
the  problem,  one  has  determined  deflections  and  residuals 
for  the  given  loading,  and  presumably  one  lias  computed 
a number  of  relaxation  patterns.  Now  the  influence  func- 
tion for  a moment  can  be  determined  as  the  deflection  of 
the  plate  due  to  a system  of  loads  grouped  near  the  point 
considered,  some  upward  and  others  downward.  A quick 
estimate  can  often  be  made  directly  of  upper  and  lower 


1 N M.  Ncwrrark.  ‘Note  :>n  Calculation  ->f  Inllurnco  Sufl.it  r»  in 
PUIm  t/>*  Ume  of  Difference  Kquaiionn.”  Journal  of  Applied  Mechanic  •, 
Vol.  &,  N«  . 2.  June  194!,  p.  A-92 

‘D.  R,  Hitrirce.  *Calc  ulu*  ing  Instrument*  »nd  Machine*,*  Univrfuly 
ot  ll):ncift  Pre»*.  Urban*.  5949,  p.  116. 


We  have  thus  ‘extrapolated*  or  ‘jumped’’  over  several 


iterations.  We  can  state  the  relations 
three  successive  steps  in  which  y - x 
Then 

also  for  the  case  of 

t It 

, z = >•  V 

u •-=  x * — X , U - X 

- 2 v * 1 X ( 3 3 1 

and 

i 

. ( * * - x ) (x  * - 2v # 

/ — — L-J  tit  1 

+ V [361 

V"'  lx  - 2 tLxj1 

/ . i i i 

The  best  value  for  the  next  iteration  is  then 


z = x t ■ t ( x - x ). 

it  ii 


(3-1 


When  we  have  only  one  value  of  i,  we  have  individual 
numbers  instead  of  a set,  and  equation  [3-]  reduces  to 
the  relation 


12 


m 


u 


z - x - (x  ' - x)2  / (x'  - 2x  ' + x)  [38] 

which  is  the  form  given  by  ilartree  tor  this  case. 

The  use  of  these  telatinnsbips  effects  a marked  im- 
provement in  convergence.  As  a matter  of  fact,  the  ex- 
trapolation leads  to  convergence  even  in  situations  ’"here 
straightforward  iteration  diverges  tapidly.  Fot  linear  sys- 
tems the  form  of  the  tesult  indicates  that  convergence  may 
always  be  expected  formally.  However  there  are  cases 
where  even  this  procedure  converges  too  slowly  to  be 
practically  useful. 

8.  Relation  of  Extrapolation  to  Relaxation 

An  interesting  relation  between  extrapolation  and  iter- 
ation can  be  derived.  If  we  use  formal  itetative  proce- 
dures in  a set  of  linear  equations,  as  in  equations  Inland 
[('].  we  find  that  we  can  relate  u and  r in  the  following 

ti  ° 

way: 

U rr.-T . / fl  l3v] 

t t tl 

Now  if  we  take  a telaxat'nn  pattern  u.  and  complete 
the  change  in  residuals  />(  from  equation  [2],  we  can  show- 
show  that  >/’.  in  equation  (35l  can  be  written  as: 

a,  = — /)./  a [40] 

t 1 n 

If  instead  of  equation  [ 33 1 we  had  chosen  to  minimize 
the  quantity 


E "b  ['•, + ' %]’• 

we  would  find  the  condition: 


In  these  relations  rr.  is  ari  arbitrary  weight  factor.  VCe 
may  take  rx . as  a ‘.  If  we  do,  we  find: 


/ «3  — 


(42) 


It  will  lie  shown  latet  that  these  equations  are  neaily 
the  same  as  those  derived  by  use  of  a slight  modifica'ion 
of  Temple’s  criterion  fot  a ‘steepest  descent*  process  :n 
relaxation. 


FORMAL  RELAXATION  PROCEDURES 

9.  General  Relations  for  Steepest  descent  Methods 

In  this  paper  consideration  is  given  only  to  formal  pro- 
cedures. Several  variations  of  steepest  descent  methods 
have  been  described  in  the  literature’,  and  the  genera! 
failings  of  these  procedures  have  also  Len  discussed  by 


’C.  Temple.  *The  Genernl  Theory  of  Kelaxalion  Mel  node  Applied  Iv 
Lin*«r  Syilrnt,"  Proc.  Koy.  S;>c.,  London,  S*r.  A,  Vol.  169.  I°!9, 
p.  '76-S00. 

H.  Couran!,  'Or,  .1  M*th*>*J  (or  th*  Solution  of  BuunJ»ty*V«]je  Problems,* 
Theodor*  von  K«rrr.*n  Anmve.'iary  Volun.e,  Cahfomih  Institute  o/  Tech- 
nology. Pa**»d<n.-»  194  1. 

D.  R.  *F.  »p.*r  iTient  al  A/ ithrnef  ic,  • Eureka.  March,  1948.  p.  13. 


Harwee  and  Coutant.  No  more  will  be  said  here  on  these 
points  which  are  still  under  study. 

If  in  a set  of  equations  such  as  [l]  we  seek  to  cortect 
or  better  the  approximate  values  x by  the  addition  of  lu 
whete  u is  a cotrection  pattern,  we  have  several  options. 
One  which  seems  reasonable  is  to  choose  n along  the 
gtadient  of  the  function 


m 


t 


--  F 
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and  to  select  t so  as  to  minimize  the  quanity  y'  rn.fr. 
vtpf.  **  ' ' 

t 

This  leads  to  the  following  process: 

(a)  For  the  given  values  of  x compute  r ftom  equations 

III. 

(b)  Operate  on  the  values  of  r as  if  they  were  e ccrr;~- 
tive  pattern,  to  determine  u.  In  this  operation,  however, 
transpose  the  coefficients  ci  the  mattix  to  use  a instead 
of  a , . That  is 

ik 


k 


[44] 


In  many  sets  of  equations  a ■»  a , or  the  matrix  is 
symmetric. 

ic)  Take  u as  a correction  pattern  and  determine  p from 
equation  [2] 

* 

E a,k  Uk  “ Pi 


(d)  Then  determine  t as  follows: 


t 


E 

E m,  p: 


[45] 


(e)  Toe  new  value  of  x is  now  x + u, 

The  procedure  described  requires  two  sets  of  opera- 
tions to  determine  the  correction  pattern,  and  it  is  there- 
fore unwie'dy. 

in  Temple’s  procedure  the  quantity  that  is  minimized  is 
the  potential  energy  of  the  system,  which  for  a symmetri- 
cal matrix  can  he  exptessed  in  the  fotrr,: 


U 


aik 


X.  X, 
t k 


E bk  xk 


[46] 


The  procedure  turns  out  to  be  as  follows: 

(a)  compute  r for  the  given  values  of  x. 

(b)  take  u.  » r ^ as  the  correction  pattern. 

(c)  compute  p * from  equation  [2]. 

(d)  compute  t fiom  the  relarion: 


t — 


147] 


A slight  modification  of  Temple’s  procedure  appears  to 
be  hictter  in  most  cases,  and  is  applicable  to  non-symmci- 
ric  ir.attices.  This  involves  arbitrarily  taking  alternatively 
two  types  of  correction  pattern. 

One  type  is 


13 


(48! 


u . « r 

t t 


and  the  other  is  a pattern  proportional  to  the  deflections, 
i.e. 

[49l 


u . - X 
l I 


In  either  case,  the  quantity  to  be  minimized  is 

I 

^ m.fr * * t p.f*  with,  in  rrc'-.t  cases,  - 1.  In  general, 
the  value  of  t is 


E w, r,  />, 
E p: 


as  in  equation  [45]) 

It  is  interesting  to  note  that  Temple’s  criterion,  equa- 
tion [47 j,  is  equivalent  to  minimizing  the  quantity 


7t2>,  + < v’-E(fr-p-) 


hoi 


10.  i'xe  o/  Orthogonality 

It  can  be  shown  that  in  general  with  the  steepest  de- 
scent procedures,  the  vectors  corresponding  to  successive 
steps  are  orthogonal,  or  nearly  orthogonal  in  successive 
pairs,  but  in  general  they  are  not  orthogonal  to  other  vec- 
tors in  preceding  steps.  It  is  entirely  possible,  however, 
to  build  up  a system  of  orthogonal  vectors  to  eliminate 
any  set  of  residuals.  .Since  mort  generality  is  obtained 
by  considering  a weighting  factor,  it  is  included  in  the 
argument. 

Suppose  we  consider  a pattern  u corresponding  to 
changes  in  residuals  p according  to  equation  [2].  A sec- 
ond pattern  f corresponds  to  changes  in  residuals  a.  It  is 
desire;!  to  find  a pattern  t corresponding  to  residuals  q 
such  that  the  following  condition  is  satisfied: 


Vie  tan  achieve  this  result  by  removing  from  q and  v 
the  portion  that  is  not  orthogonal  to  p and  u. 


Let 

Then 


;.  — v — cu 


q + cp 


and 


E m, p, - E";,  f,  E h 


[fA  1 
l > t-  J 


Hut  the  second  term  in  (52)  is  zero.  Therefore  we  have 
the  result 


E m,  p, 
E m,  p; 


and 


q = q -cp. 


r *=  t'  - cu . 


(53) 


1 54 1 


By  repeated  steps  we  can  build  up  a system  of  patterns 
and  changes  in  residuals  which  arc  mutually  orthogonal  in 
the  sense  of  equation  [51 1-  We  need  to  determine  only  as 
many  sets  of  values  as  there  are  unknowns  in  our  equa- 
tions. Then  to  solve  the  equations,  we  obtain  a set  o?  in- 
itial residuals  r corresponding  to  some  set  ol  values  of  x, 
and  proceed  to  resolve  r into  components  corresponding  to 
each  of  the  patterns.  The  coefficient  in  the  series  of  com- 
ponents corresponding  to  pn’  is  obtained  as  follows: 

Jn) 


' = t‘mP ' 


[55] 


E r,  p. 


(n) 


e p;rj  p;n) 

The  final  quantities  are  then 


E 


p.  q. 

i 1 1 t 


hi] 


6 . «=  X -r  i 6 
i ; n i 


. \rt) 


(56] 
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ABSTRACT 


A numerical  method  is  presented  for  the  solution  of 
certain  problems  concerned  with  buckling  of  flat,  rec- 
tangular plates  under  conditions  for  which  the  usual 
partial  differential  equation  of  equilibrium  can  be  re- 
duced to  an  ordinary  differential  equation. 


FUNDAMENTAL  EQUATIONS 


(5)  The  normal  stress  o > acting  in  the  x-direction  is  a 

linear  function  of  y . This  condition  is  necessary  to  in- 
sure a zero  value  of  tliroughout  the  plate. 

In  the  following  discussion  the  longitudinal  and  trans- 
verse directions  refer  to  the  x and  y directions,  respective- 
ly. 

For  the  buckling  problems  considered,  the  deflection 
surface  car.  be  represented  as  follows: 

w - Y(y)  sin  ^ Y sin  [2] 


THE.  differential  equation  governing  the  buckling  o!  a 
plate  is  given  by  Equation  (ll, 


fr'lt 

0*tv 

(fw 

P I d‘w  fr’w 

_ _2.L  4 

dx‘ 

(lx 1 fry  1 * * * V 

+ <v 

1 b 4 Pv 

N \ dx ’ * fry* 

fr’w  \ 

dx  fry) 

where  w is 

rhe  deflect 

ion,  A 

is  the  slab  stiffness,  P 

o 


parameter  proportional  to  the  loading,  />  and  (>  are  the 

* y 

compressive  buckling  forces  per  unit  of  length  acting  in 
the  x and  >•  directions,  respectively,  relative  to  the  load 
paiameter  P , and  /»  is  a similar  measure  of  shearing 

, II  IV 

force. 

A transformation  of  this  equation  into  an  ordinary  dif- 
feiential  equation  can  be  made  in  the  special  case  of  a 
rectangular  plate  with  rhe  following  edge  conditions  and 
distribution  of  loading: 

(!)  The  edges  x - 0 and  x - s arc  hinged. 

(?'  The  edges  y - 0 and  y = b may  have  any  condition 
of  support,  not  necessarily  similar. 

(3)  'ITiere  are  no  shearing  rtresses  in  the  plane  of  the 
plate  p.iraUel  to  cither  edge;  thus  r =0  everywhere  in 
the  piate. 

(4)  The  normal  stress  a acting  in  rhe  y-direction  is 
uniform  over  the  plate. 


where  Y(y)  is  a function  of  y only  a is  the  longitudinal 
length  of  one  buckle,  equal  to  x/m,  and  m is  the  number 
of  half  sine  waves  in  the  longitudinal  direction.  The 
length  a must  be  specified  before  the  solution  can  be  de- 
rived. 

'I'hen  Equation  [ 2l  is  substituted  into  Equation  [i],  and 

if  p is  set  equal  to  zero  Equation  [3I  is  obtained.  Each 
xy  1 

prime  symbol  denotes  one  differentiation  with  respect  to  >. 


The  scope  of  this  special  case  is  restricted;  however, 
many  practical  problems  can  be  simplified  to  fall  in  this 
category  and  nearly  all  available  "exact*  solutions  em- 
body the  restricted  conditions  described.  An  excellent 
discussion  of  the  classical  methods  and  a fairly  complete 
compilation  of  available  solutions  is  given  by  Timoshenko. 

This  paper  presents  a numerical  method  of  solution  of 
Equation  [3I  which  is  a combination  of  two  we  1 1-known  pro- 
cedures. Idle  Stodoln-Viaiicllo  method  of  successive  ap- 
proximations is  combined  with  rhe  numerical  procedure  of 
integration  publisher  by  Newmark.1 

Recently  Stussi'has  published  an  article  in  which  this 
same  class  of  problems  is  solved  by  a numerical  method. 
However,  there  is  only  slight  similarity  in  the  method  of 
solution  proposed  by  Stiissi  and  the  method  presented 
herein. 


1 "Theory  nf  Elastic  Stability."  by  S.  Timoshenko,  M-<Jr mw-HiI I Hook 
Cr  1 00  . . Nrw  York.  1936.. 

’ '?»on«riv*»l  ProCfilun*  lor  Cor^uting  Drflrctior.u,  Moments.  wn.l 
Hut  k!u,h  by  N.  \|.  Newmark.  Transactions  ASC.fi  t Vol.  108, 

I 043.  |>.  I 161. 

* •llrrerhnunjt  Her  iiculspannungcn  Gedrucktsr  l!echteckp*attcn."  by 
F Stusr.i.  Public. it  ion*.  Ip*  A.*'~c.  *■'*.  Bridge  s;vl  Structural  Engineer*. 

V I.  8,  W47.  |>.  237. 


EDGE  CONDITIONS 

Within  the  scope  of  the  mernod,  a number  of  different 
conditions  of  support  on  the  longitudinal  edges  can  be 
treated.  Of  these  only  the  following  three  con.inon  edge- 
conditions  are  described  in  derail: 
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(o)  Hinaed  Longitudinal  Edge 

Along  a hinged  edge  y = A.  the  deflections  and  the  mo 
ments  (hence  curvatures)  normal  :c  the  edge  are  zero. 
Therefore: 


Y **  0 and  Y"  = 0 
A A 


and 


n>  1 
n 1 h * 


/ 

— f" 
/>* 


one  obtains  Equation  (6’  from  Equation  1 4]. 


1 5b? 


(b)  Fixed  Longitudinal  Edge 

Along  a fixed  edge  y = A the  deflections  and  the  slopes 
normal  to  the  edge  are  zero.  Hence 

Y ’ « 0 and  Y ' - 0 
A A 

(c)  Free  Longitud.no!  Edge 

Along  a free  edge  the  reactions  and  the  moment*  nor- 
mal to  the  edge  are  zero.  Hence 

d'tc 

+ ft  = 0 

(9y 1 dx1 

and 

d’u/  d'w 

+ (2  - p) = 0 

dv*  dx’  dy 

where  p - Poisson’s  ratio.  When  the  deflections  defined 
by  Equation  (2]  are  substituted  into  these  last  conditi  ns, 
the  boundary  requirements  are  transformed  ir.to  the  f ,.»w- 
ing: 

h> 

Y*  771  ya 

pni/nl 

and 

h> 


(2-p)nVnt 

where  a = nh, 

GENERAL  PROCEDURE 

F.quation  [ 3l  may  be  restated,  for  convenience  in  the 
calculations,  as  follows: 

n‘  1 77*  / / n‘  I I 

y""  -2 y"-* y-p  —y" 

n * b * \ a 1 />  * * b1 


1 1 1 1 • .«#**  . :i  n 

] =Yr  4 K \o 
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The  solution  is  obtained  by  the  following  procedure: 

(a)  A set  of  deflections  Y is  assumed,  one  value  at 
each  node  point. 

(b)  With  the  assumed  deflections,  the  right  hand  sides 
of  Equations  [5a]  and  (5bl  are  evaluated.  These  equations 
are  then  solved  numerically  to  yield  the  two  component  de- 
flection functions  Y and  )'  . The  numerical  process  used 

K Is 

is  described  later.  Each  deflection  function  must  sep- 
arately satisfy  the  boundary  requirements;  then  the  sum  of 
the  two  functions  Y will  also  satisfy  the  boundary  condi- 
tions for  any  value  of  K. 

(c)  The  deflections  Y and  Y = VI  -*■  KYn  are  now  com- 

a K U 

pared.  If  they  can  be  made  identically  equal  for  a particu- 
lar value  of  K ( = P h'/X),  then  this  value  of  K defines 

o 

the  magnitude  of  the  critical  buckling  load,  and  Y*  or  Y is 

the  corresponding  configuration.  If  the  assumed  deflections 
cannot  be  made  equal  to  the  computed  deflections,  new 
values  of  / must  be  assumed  and  steps  (a)  to  (c)  must  be 

repeated  until  a desired  measure  of  agreement  is  reached. 

’J’hen  the  comparison  indicates  that  the  assumed  set  of 
deflections  is  not  correct,  a new  and  more  accurate  set  of 
detleciions  can  be  obtained  from  the  calculations.  This 
improved  deflection  configuration  can  then  be  used  as  the 
assumed  set  of  deflections  in  the  next  cycle  of  computa- 
tions. The  new  configuration  is  computed  for  the  value  of 

K which  makes  V„  t KY-  most  nearly  equal  to  Y . 

R 0 a 

The  numerical  solution  depends  upon  the  length  a. 

The  procedure  will  not  indicate  the  proper  number  of  half 
sine  waves  in  the  longitudinal  direction  to  yield  the  low- 
est critical  load  in  any  given  problem.  If  the  proper  value 
of  a is  not  obvious,  trials  with  several  values  of  m are 
necessary. 

The  functions  Y^""  and  Y^"  , expressed  by  Equations 

[5a]  and  [5b],  cannot  be  directly  evaluated  because  the 
required  values  of  Y^  ' are  not  known.  Therefore,  each  of 

these  functions  is  split  into  two  parts,  as  follows: 


14] 

in  which  K • P^bl/N  and  a ■=-  nh . The  quantity  K is  a di- 
mensionless factor  which  determines  the  magnitude  of  the  where 
given  pattern  of  buckling  forces.  The  plaie  width  b is 
divided  into  any  desired  number  of  segments  of  equal 
length  h.  Then  the  length  of  the  plate  a is  expressed  ns 
n segments  of  length  b.  It  is  not  necessary  that  n be  an 
integer. 

Using  the  notation 

77*  / 77»  f 

y""~ y"  [5a] 

M*  b‘  7i*  /j*  and 


II  I I II  1 1 


y ""  + y ' 

,rR  ,TR 


V""  -= V 

K 4 » . c 

7i*  b' 


77 3 / 

>y"' 

77 3 h' 
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T AUt  I*  OF  OPERATIONS 
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The  general  expressions  for  the  third  derivatives, 
found  by  integration,  are  given  in  Fquations  (7ai  and  [ 7bl . 


Y‘"  - ,Y‘"  * tYD"‘  - I Jy  * 2 V'  [7a] 


y"‘  ~ y * y 

D ~'rD  * 1 r‘ 


/ 

, III  _ f y.  Ill  I 

o =J‘Yo 


n*  I 


/.’  h 1 


/ 


dy-p  — y*  [7b] 

y a 


The  third  derivatives  of  the  component  deflection  func- 
tions cannot  be  completely  evaluated  unless  values  of  Y‘ 

are  known.  The  total  values  are  not  needed  except  at  a 
free  edge  or  an  edge  elastically  restrained  from  deflection. 
The  formulas  for  the  second  derivatives  are  as  follows: 


Note  that  unlike  the  third  and  fourth  derivatives,  the  sec- 
ond derivatives  may  be  directly  evaluated  from  a knowl- 
edge of  only  rhe  assumed  deflections  y . The  component 

deflection  functions,  Y^  and  Y^  , may  now  be  determined 

by  two  straightforward  integrations  of  the  total  second  de- 
rivatives,  and 

A symbolic  representation  of  the  procedure  is  given  in 
Fig.  1.  This  figure  illustrates  the  notation  used  and  gives 
a picture  of  the  plan  of  the  solution.  The  integration  sym- 
bols denote  numerical  integration.  In  step  2 of  Fig.  1 the 
quantities  ,Y^  or  , Y^""  are  evaluated.  Steps  3 and  4 

of  Fig.  1 each  indicare  one  integration  of  tiu  component 
evaluated  in  step  2.  The  total  values  of  Y"‘  are  not  de- 
termined in  Fig.  1;  if  desired  they  may  be  evaluated  in 
accordance  with  Equations  [7a]  and  [7b).  The  total  values 
of  Y"  are  computed  in  step  6.  Steps  7 and  8 are  straight- 
forward integrations.  It  can  be  seen  fr«xn  Fig.  1 that  the 
operations  used  to  determine  the  two  component  deflection 
functions  YR  and  Y are  identical. 

In  the  common  factor  column  of  Fig.  1 the  factor  for 
,y  is  c/h*,  where  c is  any  arbitrary  unit  of  deflection. 
However,  during  the  integrations  the  quantities  are  multi- 
plied several  times  by  b,  so  that  when  the  operations  are 
completed  the  factor  for  Y is  the  same  as  that  of  the  as- 
sumed deflection.  Thus  the  factors  for  all  deflections  are 
consistent.  It  is  not  necessary  to  evaluate  the  increment 
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length,  h,  in  the  course  of  the  operation. 

The  various  corrections  require*!  to  satisfy  boundary- 
conditions  are  not  indicated  in  Fig.  1,  but  are  discussed 
later. 

Consider,  for  example,  the  calculation  of  the  critical 
buckling  load  for  a square  plate,  hinged  on  four  edges,  with 
uniformly  distributed  compressive  forces  which  act  in  the 
longitudinal  direction.  The  calculations  for  an  assumed 

sinusoidal  deflection  configuration  Y are  given  in  Prob- 
lem i.  The  plate  width  b has  been  divided  into  four  seg- 
ments of  equal  length  b.  Since  the  plate  is  square,  the 
lengrh  of  the  plate  is  also  equal  to  four  increment;,  h; 
lienee  (he  value  of  n foi  this  problem  is  4.  Note  that  the 
constants  to  be  used  in  the  numerical  solution  are  given 
to  the  right  of  the  figure. 
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PROBLEM  1 


TT1 

0.  hi  7 
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n ‘ 


= 0.  3«  / 


All  Edges  Hinged 


The  presentation,  of  the  first  illustrative  picb.  em  fol- 
lows the  operations  indicated  in  Fig.  1.  Almost  every 
line  of  calculations  has  a step  number  which  corresponds 
to  the  equivalent  line  in  the  operations  table  of  Fig.  1. 
Between  steps  2 and  3 and  steps  6 and  7 are  tines  of 
computations  which  are  not  identified  with  the  operations 
table.  These  lines  contain  the  calculations  of  equivalent 
concentrations,  a basic  feature  of  the  numerical  summa- 
rion  method  as  given  by  Newmark.  A bar  above  the  sym- 
bol for  a distributed  derivative  denotes  an  equivalent 
concentrated  derivative.  Since  the  functions  lY>:"  and 
V ’ are  smooth  curves,  the  equivalent  concentrations  are 
calculated  by  the  parabolic  formulas  given  in  Fig.  2.  The 
values  of  the  odd  derivatives,  Y'"  and  Y ' , may  be  con- 
sidered as  average  values  computed  for  the  interval  be- 
tween node  points.  For  the  sake  of  brevity  the  word 
“average*  is  omitted.  Values  can  be  computed  at  the 
node  points,  if  desired. 

On  the  extreme  right  hand  side  of  the  computations  in 
Problem  1 is  a column  of  common  factors.  These  factors 
agree  with  those  in  the  operations  table  and  the  changes 
in  them  from  one  step  to  another  arc  consistent  with  the 
numerical  summation  method. 

At  the  top  of  the  calculation  is  a diagram  of  the  trans- 
verse section  of  the  plate.  This  diagram  indicates  the 
conditions  on  the  longitudinal  edges  and  defines  the  po- 
sition of  the  node  points  at  which  the  calculations  are 
being  performed. 

Only  half  of  the  plate  is  considered  rn  the  computation 
since  the  structure  and  the  deflections  are  symmetrical 
about  the  center  line.  The  derivatives  of  even  order  are 
symmetrical  about  the  center  line  and  the  derivatives  of 
odd  order  are  anti-symmetrical  about  the  center  line.  The 
boundary  conditions  considered  in  rhis  calculation  «Je 
iheretore  as  follows: 


yr"  = V'  - 0 aml  yr  -jYu^°  at  the  cd*c: 

Y'"  - Y'“  = 0 and  Y'  - Y ' * 0 at  the  center  line. 
R t J R D 


In  Problem  1 the  computations  for  steps  4 and  8 are  in- 
itiated at  the  left  edge  where  the  values  of  \Yr‘  , Y^‘ , 

Yr,  and  Y^  are  known  to  be  zero.  The  computations  for 

steps  3 and  7 are  initiated  at  the  center  line  where  , YR"'  , 

Yp"  , Yr'  , and  Y ' are  known  to  be  zero.  The  equivalent 

concentrations  given  for  the  center  node  point  are  the 
total  concentrations  of  the  function  at  :nc  center  line. 

The  equivalent  concentrations  for  the  central  node  point 
arc  tnen  split  into  two  equal  parts  for  determining  the  av- 
erage value  of  the  first  integral  in  the  segment  adjacent 
:o  the  cenual  node  point. 

The  value  of  A is  found  in  this  example  by  taking  a 
ratio  of  sums  as  indicated.  For  this  problem  the  exact 
solution  is  readily  obtained  from  Bryan’s  formula,  as 
I*  - i'A  4*  V ' (i‘ . The  value  obtained  by  the  numerical 

piocedurc  is  39.46  \/b>.  Reasonable  accuracy  can  be 
obtained  with  even  laiger  segments  in  this  problem.  For 
two  di-  isions  in  the  plate  widths  = 2)  the  critical  stress 
obtained  by  the  numerical  procedure  is  3 8.40  S/bl  and  for 
three  divisions  (n  - 3 j the  critical  stress  obtained  is 
?*>..?«  S/b\ 
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Taken  from  Fig.  5 of  Reference  (2) 


Fit-  3.  FORMULAS  FOR  EQUIVALENT  CONCENTRATIONS. 


Only  a slight  error  in  the  assumed  deflection  config- 
uration may  be  reflected  in  the  comparison  between  Y and 

Y.  Also  round-oH  errors  sometimes  make  an  exact  check 
difficult.  Fortunately,  the  value  of  K.  which  is  propor- 
tional to  the  buckling  load,  is  not  nearly  as  sensitive  as 
the  deflection  configuration  to  small  errors  in  theassume'I 
configuration.  The  required  measure  of  agreement  of  de- 
flections seems  to  be  a matter  of  judgment  since  it  de- 
pends somewhat  on  the  problem.  However,  for  a given 
number  of  segments,  the  agreement  between  Y and  Y can 

be  improved  to  any  degree  desired  by  performing  further 
cycles  of  calculations. 

The  application  of  the  procedure  to  a more  difficult 
problem  is  indicated  in  Problem  2 in  a later  section. 

ESTIMATE  OF  CRITICAL  LOAD 

In  order  to  compare  the  computed  values  of  deflection 
Y with  the  assumed  values  Y , it  is  first  necessary  to  de- 

termine  a value  of  K.  The  method  of  determining  K should 
be  such  that  a reasonable  value  is  obtained  even  when  the 
assumed  deflection  configuration  is  considerably  in  error. 
Mhen  the  assumed  deflection  configuration  is  exact,  the 
method  must  yield  the  proper  value  of  A for  the  number  of 
segments  considered.  Two  methods  are  recommended. 

The  simplest  and  most  convenient  method  of  *stimating 
K is  to  make  the  sum  of  the  errors  equal  to  zero.  The  er- 
ror at  a node  point  is  defined  to  be  the  diffctence  between 
the  computed  and  assumed  values  of  deflection.  In  this 
way  Equation  I&1  is  derived. 


K 


yd 


rsi 


A more  reliable  method  utilizes  the  principle  of  least 
squares.  In  this  method  the  value  of  K is  chosen  so  that 
the  sum  of  the  squares  of  the  errors  is  minimized.  This 
procedure  leads  to  Equation  ['/|. 


[9', 


For  most  problems  the  simple  summation  method  is 
satisfactory  and  is  preferable  because  it  is  mere  conveni- 
ent than  the  least  squares  procedure.  However,  in  prob- 
lems where  the  deflection  configuration  is  partly  positive 
and  partly  negative,  the  least  squares  procedure  is  nec- 
essary because  the  summation  method  is  not  accurate.  In 
general,  when  the  summation  method  proves  to  be  unsat- 
isfactory, the  least  squares  method  should  be  used. 

BOUNDARY  TREATMENT 

The  numerical  procedure  will  progress  in  n srr?:ght- 
forward  manner  if  the  boundary  conditions  are  such  that  a 
value  of  each  of  ,Y  , tY  " , Y ' and  Y is  prescribed  in 
advance  at  some  point  or  points.  This  is  the  case,  for 
example,  in  Problem  1 where  tY‘"  and  Y1  arc  known  to 
be  zero  at  the  axis  of  symmetry  and  tY  " and  V are  known 
to  be  zero  at  the  hinged  edge.  In  general,  other  boundary 
conditions  require  an  additional  computation. 

When  these  initial  values  are  not  all  known  in  advance, 
a boundary  value  is  assumed  where  needed  to  continuethc 
integration  procedure.  Then  at  some  later  integration  an 
inconsistent  boundary  value  will  be  encountered,  unless, 
by  chance,  the  previously  mentioned  assumption  was  cor- 
rect. !f  a boundary  value  is  inconsistent,  a correction  may 
be  made  by  finding  the  effect  on  this  boundary  value  of  a 
unit  change  in  the  earlier  assumed  value,  and  adopting  as 
much  of  this  unit  correction  as  is  necessary  to  satisfy  the 
boundary  icquirement  consistently. 

This  procedure  is  illustrated  by  Problem  2 which  pre- 
sents the  solution  of  the  problem  of  the  buckling  of  a 
rectangular  plate  subjected  to  a linearly  varying  longi- 
tudinal force.  The  two  transverse  edges  are  hinged;  one 
longitudinal  edge  is  fixed,  and  the  other  is  free.  The 
boundary  conditions  are  as  follows: 

Y1  = U and  Y - 0 at  the  fixed  edge; 

V * = ?.  75 1 h*  Y " ' and  V'  = 22.80  h'Y  " at  the  free  edge. 

The  numerical  work  is  straightforward  and  follows  FTg. 
1 and  Problem  1 except  where  explained  here.  When  val- 
ues of  tY"1  are  to  be  determined,  since  no  boundary  val- 
ues are  available,  a trial  value  is  assumed  and  the  pro- 
cedure is  continued.  Similarly  a trial  value  of  , Y muse 
be  assumed.  The  boundary  values  at  the  fixed  edge  are 
used  to  initiate  the  calculations  for  the  values  of  Trial 

Y ' and  Trial  Y. 
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PROBLEM  2 


<n  = 7.53 


0.1 755 


7 — = 0.3509 
_ > 


■n" 

0.0308 

n * 

Free  Edge  Requirements 
Y'  =3.257  h>  Y'" 

Y ' 22  80  h*  Y" 
Poisson's  Ratio  * 1/4 


O«*crl|itlon  I 5 | 

3 


o rrtfn  on 
Factors 


First  Correction  Function 


Ay,'" 

1 

1 

1.00 

1 

1.00 

1 

1.00 

I 

1.00 

» 

| 

1.00 

1 

1.00 

c/h1 

Ay," 

-5.00 

-4.00 

-3.00 

-2.00 

- 1.00 

0 

e/h' 

A?," 

-2.33 

-4.00 

O 

0 

rr\ 

1 

-2.00 

- 1.00 

-0.17 

< /h 

Ay,' 

0 

-2.33 

-6.33 

-9.33 

-11.33 

-12.33 

- 12.50 

c/h 

Ay, 

0 

-2.33 

-8.67 

- 18.00 

-29.33 

-41.67 

c 

Saco.-vl  Cor»»c*lon 

Funs  lion 

Ay," 

1.00 

1.00 

1.00 

1 00 

>.00 

1.00 

c/h 1 

Ay," 

0.50 

1.00 

1.1AJ 

1.00 

1.00 

0.50 

c/h 

Ay,' 

0 

0.50 

1.50 

2.50 

3.50 

4.50 

5.00 

c/h 

Ay, 

0 

0.50 

2.00 

4.50 

8.00 

12.50 

c 

31hen  we  denote 

• 

*,  = ""juirrd  amount  of  first  correction  function 
and  *,  = required  amount  of  second  correction  function, 
the  following  equations  express  the  free  edge  requirements. 

Y'  ♦ x,  A Y[  ♦*,  A V.'  = 3.257  />>  fV"'  ♦ x,  Ay,'"  5 
y * *,  Ay,  * x,  Ay,  = 22.80  ry"  * x,  Ay,"  ♦ x,  Ay,"; 

The  values  of  y"',  y",  y',  and  y above  are  denoted  as  "trial*  values  ir.  the  subsequent  computations.  Substituting  the  values  of 
the  correction  functions  above,  one  obtains  the  following  equations. 

c c c 

V'  - 12.  50  — x.  * 5 —x,  = 3.25  7 bl  ( Y " ♦ — x,3 
. h h h> 


Y -41. 67  c x,  * 72  50  c x,  =22.80  h>  (Y  * — x,7 

3>> 


Solving  for  x,  and  x,.  one  determines  the  following  relations. 


1 , h 

74.13  x,  - — (Y  - 22.80  h'  Y ) * 2060— (Y‘  - 3.257  h 1 y'") 
c c 

7 „ b 

23.52  x,  =— fy  - 22.80  y 7 - 2.644  ~(Y'  - 3.257  .■>'  Y'") 
c c 
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Common 
F artor* 


VTfi"1Y* 


Calculation  zi  Yn 
I K . 


at  free  end  Y * 5.7  <*//>  (calculated  from  previous  cycle). 


■V" 

-0.052  -0.172 

-0.326 

-0.487 

-0.653 

c/h 4 

rv'"' 

-0.058  -0.1. ’5 

-0.327 

-0.487 

-0.299 

c/h' 

Tr  i al , VR  " ' 

- 1.153  -1.212 

-1.387 

-2.714 

-2.201  -2.500 

c/h ' 

Trial  ,YR"  0.069 

-1.086  -2.298 

- 3.68' 

-5.399 

-7.600 

c/h 

0.596  1.965 

3. '20 

5.544 

7.439 

c/h' 

Trial  YR"  0.0(>8 

-0.490  -0.333 

0.033 

U.  145 

-0. 161 

c/h 

Trial  YR"  -0.089 

-0.430  -0.315 

0.014 

0.110 

-0.012 

c/h 

Trial  V'  0 

-0.089  - 0.519 

-0.834 

-0.820 

-o. 

710  -0.722 

c/h 

Trial  VR  0 

-0.089  -0.6>08 

-1.442 

-2.262 

- 2.972 

c 

at 

the  free  edge 

YH 

-22.30  h 1 VR"  - -2.972  <-  * f 22.  80  ) 60.  .'6/  < ) 

= 0. 699  r 

it*  / r / c \ <* 

Vo  " = .Vo"'  * 2 v'*=  -2.500  — * (0.  5509)  (5.7  ) - -0.500  

H K n'  h'  a h ■ \ h' ) h' 

VR  - 3.257  h'  V "’  - -0.  722  — ‘ 6 3.257)  (o.  500  — ] = 0.  906  — 

6 \ 6 / 6 


Therefore 

*,=  0.0346  and  *,=-0.072/ 

0 -0.21  -1.05 


*0.0  346  Ay, 
-0.072/ Ay, 


at  free  end  VR  = Trill  VR  * 0.0 3«6  Ay,'  - 0.072/  Ay,'  = -1.52  c/h 


Calculation  of 


• p* 

0 

o._*o 

0.40 

0.60 

0.80 

1.00 

Y " " «=  Y,;" 
| 1 rP  Y0 

0 

0.06 

0.39 

1.12 

2.22 

3. '2 

c/ 

- »/. 
rD 

0 08 

0.42 

1.15 

2.25 

1.59 

c/h 

Trial  Vp'" 

-2.24 

-2.16 

- 1.74 

-0.59 

1.66 

3.25 

c/h 

Trial  Vp" 

6.77 

4.53 

2.37 

0.63 

0.04 

1.70 

c/h 

Trial  Vp" 

3.01 

4.54 

2.40 

0.'3 

0.23 

0.48 

c/h 

Trial  Vp 

0 

3.01 

7.55 

9.95 

10.68 

10.91 

11.39 

c/h 

1 ri  al  Vp 

0 

3.01 

10.56 

20.51 

31.19 

42. 10 

c 

At  tne  free  edge 

VD  - 22. SO  S’  VR"  = 42.  /0  c - 622.80)  6/.  70  cj  = 3.  34  c 

V0'  - 3.237  h ' V0'"  = //.  39  — - (3.  257)  (3.25  — ] = 0.80  — 

/>  \ h f b 
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r'RODLEM  2 (continued) 


(P)  = 0.  66  3 = 16.6 

° rf  h>  b1 


The  boundary  requirements  at  the  free  edge  are  not  in- 
itially satisfied;  therefore  the  assumed  values  of  ,Y"' 
and  ,Y”  are  not  correct.  Note  that  at  the  free  edge  the 
boundary  value  of  Y ' is  a function  of  Y '“  . However, 
only  " is  given  in  the  tabulated  computations.  In  this 

case  Y "'  is  computed  from  Equation  7(a),  which  in  turn 

requires  the  slope  of  the  assumed  deflection  curve  at  the 
free  edge  to  be  known.  Therefore,  when  one  of  the  bound- 
aries is  a free  edge  the  slope  Y1  at  that  edge  must  also 

be  assumed  and  the  final  computed  slope  must  agree  with 
this  value,  just  as  the  assumed  and  computed  deflections 
must  agree  in  order  for  the  solution  to  be  valid.  I:  is  us- 
ually more  accurate  to  take  the  slope  of  the  curve  at  the 
boundary  from  the  preceding  computation  rather  than  to 
use  an  interpolation  formula.  For  this  ease  the  difference 
is  slight  since  the  slope  Y1  found  by  the  use  of  either  a 

three  point  or  four  point  approximating  polynomial  is 
5.5  c/h  as  compared  with  the  value  of  5. 7 c/h  found  from 
the  previous  cycle.  The  latter  is  more  accurate. 

Two  correction  configurations  are  required,  one  for  each 
assumption.  Each  correction  configuration  influences  both 
boundary  requirements.  The  correct  proportions  of  the  first 
and  second  correction  configurations  (Ay,  and  Ay,  , re- 
spectively) may  he  determined  from  the  equations  presented 
just  Oelow  the  computations  for  the  configurations. 

The  values  of  the  required  corrections  are  added  to  the 
values  of  Trial  Y to  yield  the  correct  deflections;  similarly 
the  correct  end  slopes  are  found.  Both  and  YD  axe  cor- 
rected independently  of  each  other. 

The  value  of  K is  found  to  be  0. 6o),  which  gives  a val- 
ue of  the  critical  buckling  stress  (at  the  free  edge)  of 
/ 6.6  N/b 1 . 


*S«.-  -Bounds  and  Convergence  of  ReU*«t»on  and  Iteration  Proce- 

lure#.*  by  N.  M.  Newnmrk.  the**  Proceeding*. 


CONVERGENCE 

In  most  problems  the  revolting  deflection  configuration 
y is  a better  approximation  to  the  true  buckling  configuia- 
tion  than  the  assumed  deflection  shape  Y . In  this  case 

the  process  is  said  to  converge  toward  the  true  solution. 

On  the  other  hand,  in  limited  classes  of  problems  Y may 
be  further  than  Y from  the  true  solution,  and  in  this  case 

the  process  is  divergent.  A problem  in  which  the  piocess 
is  convergent  may  be  handled  most  simpiy  by  using  the  re- 
sultant deflection  configuration  Y of  the  last  cycle  for  the 
assumed  configuration  Y^  of  the  next  cycle.  This  obvious 

procedure  will  not  be  efficient  when  the  process  converges 
slowly;  m this  case  extrapolation  as  described  below  is 
desirable.  VChen  the  process  diverges,  an  extrapolation 
may  be  used  to  obtain  an  answer  in  many  cases.  In  others, 
no  solution  is  possible  with  successive  approximation 
techniques  of  this  type. 

It  is  not  a simple  matter  to  determine  in  advance  whether 
or  not  the  successive  approxim  itions  converge.  However, 
from  experience  in  solving  a variety  of  problems  it  maybe 
stated  qualitatively  that  if  there  are  tensile  forces  in  the 
plate,  the  larger  the  amount  of  tensile  force  the  greater 
will  be  the  tendency  to  diverge.  It  may  also  be  stated 
that  if  divergence  occurs,  the  successive  solutions  will 
generally  oscillate  about  the  true  solution.  These  fea- 
tures are  similar  to  those  which  apply  in  other  types  oi 
buckling  problems  in  which  successive  approximation 
techniques  are  used. 

The  recommended  extrapolation  procedure*  is  applied 
as  follows:  Let  successive  deflections  at  any  point  be 
designated  as  A.  B,  and  C,  respectively,  where  A is  the 
first  assumption,  B follows  from  A,  and  C from  B.  Then, 
the  extrapolated  value  A is: 

AC  - 

T 

A - C - 211 
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If  this  formula  is  used,  it  must  be  applied  separately 
to  the  deflections  at  each  point.  Thus,  after  two  steps 
in  the  calculation,  which  yield  values  of  H and  ?t  all 
poinrs  corresponding  to  assumed  values  A,  one  computes 
the  values  of  A,  and  repeats  the  calculation.  Ordinarily 
one  may  have  to  go  through  this  operation  several  times 
in  problems  which  diverge  rapidly. 

Problems  in  which  divergence  occurs  present  a diffi- 
culty which  becomes  ‘mmediately  obvious  to  the  computer. 
He  may  have  to  Use  other  techniques  to  solve  the  problem 
if  the  successive  approximations  procedure  does  not  lead 
to  an  answer  either  directly  or  with  the  extrapolation  pro- 
cedure described. 

SUMMARY 

The  procedure  described  is  useful  for  a wide  variety  of 
loadings.  With  some  modifications  the  problem  of  the 
buckling  of  ribbed  plates  is  easily  solved.  In  addition, 
many  other  practical  problems  which  do  not  satisfy  all  of 
the  requirements  of  this  method  can  be  treated  by  either 
simplifying  the  problem  or  by  waiving  some  of  the  limita- 


tions of  the  method  which  may  not  be  important  for  a par- 
ticular case.  Extension  of  the  procedure  to  problems  which 
do  not  satisfy  the  limitations  of  the  method  requires  en- 
gineering judgment. 

The  results  are  approximate  but  can  be  made  as  ac- 
curate as  is  desired  by  taking  a sufficient  number  of  di- 
visions in  the  plate  width.  *s  the  number  of  divisions  is 
increased  the  accuracy  of  the  solution  increases  and  also 
the  amount  of  work  required  to  effect  a solution  increases. 
The  question  of  how  fine  a network  to  use  depends  some- 
what on  the  problem  and  can  best  be  decided  by  experi- 
ence. Ordinarily,  however,  for  most  practical  problems  nc 
more  than  five  or  six  divisions  in  the  plate  width  are 
needed. 
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INTRODUCTION 

A STRESS  solution  of  most  radially  symmetric  problems 
is  difficult  and  tedious  by  formal  mathematical  pro- 
cedures. In  this  work  a numerical  method  is  developed 
that  can  be  readily  applied  to  solve  such  problems.  The 
particular  problem  selected  to  illustrate  the  method  is 
that  of  the  restrained  cylinder  under  compression.  The 
method  requires  a framework  of  pin-connected  bars  which, 
for  similar  load  conditions,  will  approximate  the  deforma- 
tion of  the  solid  body.  This  method  of  solving  stress 
problems  is  referred  ro  as  the  lattice  analogy  me I bod.  The 
arrangement  of  the  bars  in  any  element  of  the  framework 
is  arbitrary.  Normally  a pattern  that  is  simple  and  suit- 
able for  the  type  of  problems  being  considered  is  adopted. 

It  is  apparent  that  the  smaller  the  size  of  rite  mesh  of  the 
framework  the  better  the  approximation  of  the  deformation 
of  the  solid  body  by  rhe  framework. 

To  determine  the  deformation  of  the  framework  it  is 
necessary  to  know  the  area  of  each  bar.  These  areas  are 
derived  from  the  requirement  that  the  distortion  of  an  ele- 
ment of  the  solid  body  and  of  the  framework  be  the  same 
for  similar  load  conditions.  Once  the  bar  areas  are  es- 
tablished the  relationship  between  the  force  and  th<"  end 
displacements  of  each  bar  can  be  formulated.  Then  from 
the  equilibrium  conditions  of  the  joints,  displacement 
equations  are  developed.  In  these  eqjations  a component 
displacement  of  a central  joint  is  related  to  the  two  com- 
ponent movements  of  tlf  immediate  surrounding  joints. 

The  set  of  simultaneous  equations  that  arise  in  a partic- 
ular problem  is  solved  either  by  an  iteration  or  relaxation 
technique.  Fc-  the  problems  solved  in  this  work  an  iter- 
ation process  was  used. 

The  advantages  gained  in  the  use  of  a lattice  analogy 
method  are  often  overlooked.  Once  the  size  A.  the  mesh 
of  the  framework  has  been  selected  the  only  approxima- 
tion involved  in  a solution  of  a problem  is  the  degree  of 
precision  to  which  the  numerical  work  is  carried.  U’ith 
the  use  of  finite  differences  besides  the  approximations 
due  to  the  si.x-  of  the  mesh  and  the  precision  of  the  nu- 

•Numbers  tr  pa/enlhrtc*  r»fer  10  referer.ee*  li*te<1  in  th?  bibliography 
at  the  rr.d  of  the  paper. 


merical  work,  there  is  the  approximation  of  the  differen- 
tial equation  by  differences.  Also,  the  difficulties  that 
arise  at  boundaries  with  a finite  difference  procedure  for 
example,  at  a fight  a.ig.e  or  re-entrant  corner  are  not  pres- 
ent in  the  lattice  analogy  procedures.  The  framework  also 
lends  itself  to  the  solution  of  elastic-plastic  problems. 

As  with  a finite  difference  procedure  where  the  solu- 
tion of  a problem  can  be  set  forth  either  in  terms  cf  a 
stress  function  or  i/i  terms  of  displacements,  the  lattice 
analogy  method  can  be  formulated  in  rerms  of  bar  forces 
(1),  (2)*  or  in  terms  of  displacements  (3),  (4).  The  meth- 
od developed  in  this  paper  deals  with  displacements.  The 
main  advantages  in  using  displacements  are  (1)  the  bound- 
ary conditions  can  be  prescribed  either  in  terms  of  dis- 
placements or  forces  (2)  the  body  can  be  eithei  singly  or 
multi-connected.  The  numerical  displacement  procedure 
of  analysis  appears  to  offer  the  only  means  of  solving 
stress  problems  for  which  the  body  is  multi-connected  and 
the  boundary  conditions  are  prescribed  forces  and  dis- 
placements. 

The  assumptions  made  in  the  classical  rheory  of  linear 
elasticity  apply  to  this  work. 

NOMENCLATURE 

The  following  nomenclature  is  used  in  this  paper: 

r,  0,  z - cylindrical  coordinates 

a . a.,  o - normal  stresses  in  the  r,  0,  z directions 

r ~ Z 

r - sheating  stress  in  rhe  r-z  plane 
( , t. . f = normal  strains  in  the  r,  0,  z directions 

T Z 

v - detrusion  in  the  r-z  plane 

’’7. 

i'(ii  *(  = dilatation 
^ z 

F ■=  modulus  of  elasticity 
C - shear  modulus 
p - Poisson’s  ratio 

u,  iv  = radial  and  axial  displacement  components 
A - area  of  a bar 
e = strati,  in  a I ar 
s - stress  in  a bar 
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/•'  - force  in  ! .a 
A = length  parameter 
n -r  :..;cge r 

R,  7 ^ normal  force  concentrations  at  a joint 

V - shearing  force  concentrations  at  a joint 

a,  a,  A - arbitrary  constants  used  with  expressions 
° for  the  bar  areas 

r,  f , r.7,  arbitrary  constants  used  with  the  ‘second- 
ly, /,  p ary"  force  systems 

u =^—  fa  - o J*  - f oa  - a 7 - ( a -«)’  i ' I 1 

' 2 \ ' ° z z t rrl 

- intensity  of  stress 

THEORY 

The  stress  solution  for  a radially  symmetric  body  in 
obtained  most  convenientl)  with  reference  to  a system  of 
cylindrical  coordinates.  A [>oint  P (r.  0,  z)  is  defined  by 
• rs  radial,  tangential  and  axial  coordinate  icspectively. 

I'he  radial  and  axial  displacements  of  the  distorted  body 
are  designated  by  u and  w,  respectively.  The  normal 
stress  components  designated  by  a,  oc,  and  o act  in  the 

radial,  tangential  and  axial  direction  and  the  shear  stress 
.-  acts  in  the  r-z  plane.  The  remaining  stress  components, 

due  to  symmetry,  are  zero.  Based  on  the  assumption  of 
small  displacements  the  normal  strain  components  and  do- 
trusion  are  respectively 
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that  the  differential  equaHr>"s  to  be  satisfied,  subject  to 
boundary  conditions,  by  the  displacement  components  are 
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However,  with  either  the  differential  equations,  Equation 
[3],  or  those  in  terms  of  stress  components  a solution  of 
a particular  problem  is  intractable  by  formal  mathematical 
procedures  n.-.d  zr,  app.-oximaie  numerical  method  of  an- 
alysis has  to  be  resorted  to.  Such  a method  is  developed 
in  the  follow  ing  work.  Since  it  involves  displacements  it 
can  be  considered  as  the  analog  of  the  integration  pro- 
cedure that  is  necessary  to  integrate  equations  [3]. 

A suitable  pattern  of  hars  is  shown  in  Fig.  1.  In  an 
axial  plane  the  framework  is  built-up  cf  square  lattices 
with  crossed-diagonals.  A cross-section  is  composed  only 
of  radial  and  tangential  bars.  To  establish  the  har  areas, 
the  ht-ess  systems  required  to  produce  identical  deforma- 
tions of  the  framework  and  solid  body  are  considered. 

Since  the  deformation  cf  an  element  can  be  described  fun- 
danentally  in  terms  of  three  strain  components,  the  de- 
formations considered  are  those  associated  with  a uniform 
strain  in  the  radial  and  axial  direction  and  a detrusion. 
Thus  for  the  solid  body  a uniform  strain  in  the  axial  di- 
rection of 

r -=  / 


For  an  elastic  material  which  obeys  Hooke’s  Law  the 
stress  strain  relations  arc  given  bv 
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The  solution  of  a particular  problem  requires  either  a 
determination  of  the  stress  components  or  the  displace- 
ment components.  For  example,  it  can  be  readily  shown 
from  a consideration  of  (1),  and  (2)  and  the  equations  of 
equilibrium 


F.«.  1.  FRAMEWORK  OF  UARS. 
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is  Droduted  by  the  stress  disiribution.s 


n r 


, Uc  

<1  * fi)  (!  ~ 2,0 

(1  -fi)  F. 
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r (1-10(1-2;.'.) 


r = 0 

r 7 


A uni  for  .ft  strain  in  the  radial  direction  that  is 
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requires  the  stress  distribution 
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n * ;o  n - ?)0 
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For  a dot  fusion 
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where  A , a and  a are  arbitrary  constants  A and  <p  are  de- 
fined in  Fig.  1 and  n is  an  integer.  All  the  areas  with  the 
exception  of  those  of  the  tangential  bars  are  dependent  on 
{ 4 1 the  radial  distance. 

Other  selections  fot  the  bar  areas  are  possible,  how- 
ever, with  this  selection  the  resulting  bar  areas  and  dis- 
placement equations  to  be  developed,  take  on  simple 
expressions  convenient  for  numerical  computation. 

Because  of  the  selected  bar  pattern  and  areas  it  was 
found  necessary  to  introduce  into  the  fcrce  system  for 
each  type  of  deformation  a ’secondorv"  system  of  bar 
forces.  Tithput  the  use  of  these  ‘secondary"  systems  the 

coefficients  A , a and  a in  Equations  {?]  would  turn  out  to 
o 

be  functions  of  radial  distance,  thus  complicating  the  form 
of  the  displacement  equations.  Further,  the  use  ot  such 
systems  with  a square  lattice  having  crossed-diagonals 
permitted  a determination  of  the  bar  areas  in  terms  of  an 
arbitrary  value  of  Poisson’s  ratio.  In  previous  uork,  [4], 
secondary  bar  patterns  had  to  be  introduced  when  Poisson’s 
ratio  differed  from  1/3.  The  use  of  secondary  bar  patterns 
leads  to  complications  and  is  to  be  avoided. 

Without  loss  of  generality  and  to  avoid  carrying  unt-s- 
[^1  sential  constants  A and  o are  taken  as  unity  hereafter. 

Considering  the  framework,  the  forces  due  to  a uniform 
axial  strain  are  shown  in  Fip.  2.  Component  forces  maA 

o 

and  qmaA  are  the  resultant  radial  components  of  the  forces 

in  the  tangential  bars.  The  forces  in  the  radial  and  tan- 
gential bars  given  by  k(n  + 1/2)  A an d are  due  to 

the  ’secondary"  force  system  associated  with  a uniform 
axial  strain.  From  the  statics  of  joints  n and  e,  Fig.  2, 
the  equivalent  uniform  stresses  are 


the  stresses 


are  fc^uiied. 

Before  considering  states  of  deformation  similar  to  the 
foregoing  for  the  framework,  it  is  expedient  to  choose  the 
form  of  the  expression  for  the  area  of  each  bar.  It  is  as- 
sumed that  the  area  of  each  bar  is  proportional  to  'he 
corresponding  cross-sectional  area  of  a solid  element  with 
height  /v,  radial  length  A and  average  width  < n - 12)  A6. 
For  example,  the  area  of  the  axial  bar  is  proportional  to 
the  axial  cross-sectional  area  tj A*  of  the  solid  element, 
and  likewise  for  the  remaining  bars.  Hence  the  areas  for 
the  radial,  axial,  tangential  and  diagonal  bars  are 


A - 

r 


A = mV  cS  A 

7 O 


A.  = a .V  /t_ 


r.,;.  j.  oar  forces  due  to  cm:  orw  axiai.  strain- 
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c 


From  the  condition  of  radial  equilibrium  of  join: 


a 

rna  k * — (9l 

2 


In  a similar  manner  the  bar  forces  corresponding  to  a uni- 
form radial  strain  are  shown  in  Fig.  3.  In  this  case  the 
"secondary*  system  is  made  up  of  the  force  cnA  in  the 
a 0 

axial  bars  and — (n  + l/2)tA  in  the  diagonal  bars.  A 
4 0 

summation  of  forces  at  joints  n arid  e yields 


s ■=  lev  — (/+/,)  A 

z L 2 J 0 

* 

Hi  *■  — + qa  4 — * — j ( I -*  rjij  ( 10] 


Fig.  J.  OAR  FORCES  DUE  TO  UNIFORM  RADIAL  STRAIN. 

At  joint  r the  sum  of  the  radial  forces  gives 


• 


The  force  * a f>A  is  to  be  associated  with  the  *sec- 

~ 2 ° 
ondary*  force  system. 

In  the  above  work  the  constants  k,  m,  q,  c,  t,  and  p are 
to  be  determined  along  witl.  a,  a and  A . These  constants 

are  evaluated  by  equating  corresponding  equations  of  stress 
for  similar  states  of  deformation  in  the  solid  body  and 
framework.  Hence  from  Equations  [4]  and  [bl,  (5l  and  [lOl, 
[/>!  and  [12]  and  simplifying  by  means  of  Equations  l9l  and 
(ill  and  taking  q = l 2,  one  obtains  the  following  expres- 
sions for  the  constants: 


2(1  - (0 

( I * p)  ( l - 2p)  (2  r a) 

2/i  < 1 - at  - a 
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2(l-i t) 
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Corresponding  to  a detrusion  of  — ■*  — in  the  solid  body 

2 T 

the  strain  in  the  diagonal  bars  is 

1 

Ti 

and  all  other  bar  strains  are  zero.  The  bar  forces  for  this 
case  are  shown  in  Fig.  4.  The  equilibrium  of  joints  n and 
e requires  that 

9 n = S (rt  - Ij  - a ( l + t>)  A 112! 

rz  zr  o 


a = / - k - n (l3l 

a 

n,a  « k - — 

2 

c-k 

l - k - a 

P - 

a 

Since  a total  of  nine  constants  are  to  be  determined 
from  the  above  seven  expressions  two  of  them  can  he 
chosen  arbitrarily.  For  th-*  work  to  follow  a value  of  1/4 
is  taken  for  Poisson’s  ratio.  For  simplicity  of  results  a 
was  taken  as  one.  Accordingly,  the  values  for  the  con- 
stants are  k.  - p = c - n,  a - 2,  l = I,  m - 1/4.  4 - 4/3 
and  as  previously  assigned,  q =■  1/2. 
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^itn  these  results  the  bar  areas  become 
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It  is  interesting  to  note  that  for  plane  strain  problems 
(fi  = f/41  the  bar  areas  of  a framework  of  square  lattices 

4 4 

with  crossed-diagonals  would  be  — and  for  the  hori- 

5 

zontal  and  vertical  bars  and  the  diagonal  bars  respectively, 
For  convenience  and  summary  the  bar  forces  and  stresses 
for  each  type  of  deformation  considered  in  the  develop- 
ment of  the  bar  areas  are  given  in  Figs.  5,  6 and  7. 


Fi*.  7.  DFTRUSION  r 


J/r. 


wotk  are  fractional  parts  of  the  corresponding  interior  bar 
arias.  Thus  the  radial  and  tangential  bars  in  a radial 
plane  have  one-half  the  area  of  the  same  interior  bar.  At 
a coinei  the  tangential  bar  area  is  one-quarter  of  the  same 
interior  bar.  At  an  axiai  boundary  a distance  n A from  tiic 

7 

origin  the  axial  bar  has  an  area  of — f n - i/2).  This  is 

5 

consistent  with  the  area  4/5  n of  the  interior  bar  at  a dis- 
tance nA  from  the  origin  since  it  is  considered  that  the 

area  4/5  n is  made  up  of  two  bars  with  areas 
arid  — (n 


'*iH) 


1 2).  From  similar  reasoning,  the  area  of  each 


axial  bar  at  nA  - 0 is  1/5. 

Having  established  values  for  the  bar  areas  and  stress- 
strain  relations,  shown  diagrammatical  ly  in  Figs.  5,  6 and 
7,  one  can  proceed  to  a derivation  of  the  displacement 
equations.  Tor  a general  interior  joint  as  shown  in  Fig.  8 
forces,  stresses,  etc.,  of  all  bats  radiating  from  joint  r 
are  denoted  by  the  It-. ter  or  letters  of  the  joint  remote  from 
r.  For  example  s stands  for  the  stress  in  bar  r-c,  or  /• 

r nr 

is  the  force  in  har  c-ne.  However,  the  symbols  t.  anu  u 


Fl*.  5.  UNIFORM  AXIAL  STRAIN,  e 


Fig.  6.  UNIFORM  RADIAL  STRAIN,  e = I. 

From  Fig.  5 it  can  be  seen  that  each  unit  of  axial 
strain  requires,  besides  the  expected  bar  stresses,  a 
"secondary*  stress  of  1/4  in  the  tangential  bars  to  main- 
tain radial  equilibrium.  Likewise,  Fig.  6 shows  a "sec- 
ondary" stress. of  1/2  is  required  in  the  diagonal  bars.  In 
developing  the  displacement  equations  these  "secondary" 
stress  systems  are  taken  into  account. 

The  area  of  the  bars  at  the  boundaries  of  the  frame - 
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M-J*- 


refer  to  displacement s ar  a joint  so  that  z implies  the 

radial  dis  >laccment  of  joint  sc.  Consideting  then  a de- 
formed fra  me  wo  lk  and  using  the  stress  strain  relations 
derived  one  can  write  for  the  bar  stresses 
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The  equilibrium  of  joint  c,  Fig.  7 requires 


In  a manner  similar  to  the  above,  displacement  equa- 
tions can  Le  derived  for  the  various  boundary  joints  that 
occur  in  a problem. 

To  solve  a particular  problem  one  must  first  select  the 
size  of  the  mesh  of  the  framework.  Also,  for  that  portion 
of  the  surface  of  the  solid  body  subjected  to  stress  a sys- 
tem of  equivalent  force  concentrations  must  be  placed  on 
the  boundary  of  the  framework.  The  linear  set  of  displ ace- 
115]  ment  equations  that  occur  in  a particular  problem  car.  best 
be  soived  by  either  an  iteration  or  relaxation  procedure. 

A discussion  of  such  procedures  or  the  details  as  to  how 
a solution  of  a problem  is  carried  out  r.eed  not  be  given 
here. 

From  the  determination  of  the  joint  displacements  the 
force  in  each  bar  can  be  established  by  multiplying  the 
appropriate  stress  equation,  Equation  [151,  by  its  corre- 
sponding bar  area.  Equation  [14].  The  average  uniform 
stress  components  at  each  joint  are  then,  obtained  from  a 
cons' deration  of  sections  as  shown  in  Figs.  9 and  10. 

That  is  for  an  axial  section  the  average  radial  stress  at 
joint  c is 
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where  ar.-J  V are  the  component  body  forces  acting  at 

c.  Multiplying  the  expressions  in  F.quations  [lj]  by  their 
appropriate  area  as  given  by  Equation  (l4l  and  substitut- 
ing into  Equations  [ 16]  and  simplifying,  one  obtains  for 
the  radial  and  axial  displacement  equations  of  the  joint 
c, 
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The  tangential  stress  is 


F 

rt 
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From  a radial  section,  the  average  axial  stress  is 
given  by 
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Fig.  10.  SECTION  R-R  FOR  EVALUATION  OF  AXIAL  STRESS. 

7. 

c 

=. — 

n 

To  determine  the  shearing  stress  at  joint  c the  horizontal 
forces  above  or  below  section  R-R  in  Fig.  9,  or  the  verti- 
cal forces  to  left  or  right  of  section  7.-7.  in  Fig.  10  are 
considered.  Thus  for  joinr  c 
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(or  section  R-R,  or  equally  well  from  section  7-7, 
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When  F and  A are  raken  different  from  unity  it  can  be 
readily  shown  that  all  stress  components  have  to  be  mul- 
tiplied by  the  factor  F./X. 

The  lattice  analogy  method  of  solution  developed  in 
the  foregoing  work  was  applied  to  the  problem  of  the  end 
restrained  cylinder  subjected  to  compression.  The  re- 
sults of  the  analysis  are  discussed  in  detail  and  com- 
pared with  existing  solutions. 

END  HESTRAINED  CYLINDER  UNDER  COMPRESSION 

The  problem  of  the  cylinder  subjected  to  compression 
was  first  solved  in  a partial  manner  by  Filon,  (5),  in  1902 
and  later  by  Pickett,  ((>),  in  1944.  The  solutions  in  either 
case  were  not  complete.  The  solution  given  b)  Filon  did 
not  sacisfy  all  of  the  boundary  conditions.  Pickett  was 
unable  to  secure  without  difficulty  the  stress  distribution 
near  the  outer  contact  edges  of  the  cylinder  ard  rigid 
blocks.  Except  for  the  axial  stress  a comparison  of  these 


two  solutions  showed  large  differences  between  the  val- 
ues of  the  stress  components.  Because  of  its  engineer- 
ing significance  it  was  felt  that  a study  to  c.neck  the 
foregoing  investigations  and  to  determine  the  maximum 
values  of  the  stress  components  at  the  outer  edges  of 
contact  was  necessary. 

Consider  a cylinder  compressed  between  rigid  rough 
end  blocks  as  shown  in  Fig.  11.  The  modulus  of  elas- 
ticity of  the  material  is  taken  as  one  nnd  Poisson's  ratio 
is  1/4.  The  boundary  conditions  are: 


r 


Fit.  11.  END  RESTRAINED  CYLINDER. 

1.  No  stress  across  the  curved  surface,  that  is, 

o=0  at  r=J 

r 

r » 0 at  raj 

rz 

2.  On  the  end  planes  at  z = J, 

u » 0 for  all  values  of  r 

3 . u-  - J i at  z = + } 

Condition  1 does  not  hold  at  the  edges  z = * J,  r = J.  At 
these  edges  the  radial  and  shear  stress  arc  discontinuous. 

For  a framewoik  with  a rncsn  parameter  of  A « 1/2  the 
component  displacements  are  given  in  Fig.  12.  The  stress 
distributions  computed  from  these  results  along  with  a 
plot  of  lines  of  constant  stress  are  shown  in  Figs.  13.  4, 
15,  and  16.  A comparison  of  the  stresses  obtained  by  the 
lattice  analogy  method,  wirh  those  of  Pickett  and  Filon  is 
illustrated  in  Fig.  17.  It  can  be  seen  that  the  lattice  an- 
alogy solution  is  in  good  agreement  with  that  due  to  Pick- 
ett. Because  of  difficulties  in  rhe  convergence  of  the 
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Klc.  13.  DISPLACKMENTS  FOR  END-RF.STR  At  NF.D  CYLINDER 

B » 1.  ,*  = 1'i,  \=  1/2.  AVERAGE  COMPRESSIVE 

STRr.SS  - n - - 0.  U8. 

« 1 


Fi*.  13.  RATIO  OF  AXIAL  STRESS  o TO  AVERAGE  COMPRESSIVE 

z 

STRESSES  cr  * - 0.  \48. 
o 

series  solution  for  the  stresses,  Pickett  was  unable  to 
establish  maximum  values  for  the  stresses  at  the  outer 
edges  of  contact.  The  maximum  ratios  of  axial,  radial, 
tangential  and  shear  stress  to  average  compressive  stress 
determined  by  the  lattice  analogy  method  are  1.70,  0.64, 


4. 


Flf.  16.  RATIO  OF  SHEARING  STRESS  t to  AVERAGE  COk»>RES- 

SIVE  STRESS  cr  ="0.348.  rX 
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im<4 

film 

ruft* 

UttlM  twKo  IWM  


F If.  17.  COMPARISON  OF  STRESS  COMPONENT  RATIOS  ALONG 
z = 3.  AS  DETERMINED  BY  FILON,  PICKETT  AND  BY  THE  LATTICE 
ANALOGY  METHOD 


la*  ad 

Mam0  M IlMl 


ath  M QmH 


F‘«.  IS.  COMP  PRISON  OF  STRESS  COMPONENT  RATIOS  ALONG 
X * 3 FOR  CYLINDER  COMPRESSED  BETWEEN  ROUGH  AND  SMOOTH 
RIGID  END  BLOCKS.  COMPUTED  BY  THE  LATTICE  ANALOGY 
METHOD. 

that  the  maximum  values  would  differ  little  from  those 
given. 

Since  the  solution  given  by  Filon  satisfied  the  bound- 
ary condition  2 only  at  r = 3,  the  shear  distribution  shown 
in  Fig.  17  was  required.  This  difference  in  the  boundary 
conditions  had  relatively  little  effect  on  the  magnitude  of 
the  axial  stress,  but  it  materially  changed  the  radial  and 
tangential  stresses. 

Tc  show  the  effect  of  different  distributions  of  shear 
stress  over  the  end  faces,  a cylinder  constrained  by  rigid 
blocks  with  smooth  surfaces  was  considered.  The  bound- 
ary conditions  are  the  same  as  those  ot  the  previous  prob- 
lem except  for  the  end  planes,  condition  2 is  u = 0 at 
r * 3 only,  and  r = 0 for  all  the  vaiues  of  r.  A compari- 
son of  the  stress  components,  at  z * 3,  with  those  of  the 
previous  solution  is  shown  in  Fig.  18. 

The  average  ..ompressive  icress  for  the  cylinder  com- 
pressed between  rough  rigid  blocks  was  0.348.  Based  on 
this  value  the  modulus  of  elasticiry  as  in  a standard com- 
pression test  would  be  given  is  1.05.  Hence  to  obtain  the 
true  modulus  it  is  necessary  to  correct  that  obtained  from 
a standard  compression  test  by  the  factor  0.953.  Using 
the  Prager-Synge  method,  Edelman,  (7),  obtained  a correc- 
tion factor  of  0.948. 

Shown  in  Fig.  19  is  a plct  of  the  ratio  of  the  intensity 
of  stress  to  the  shear  mou'ilu:..  It  is  interesting  to  note 
that  after  initial  plastic  ilow  in  a localize.!  region  near 
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the  outer  edge  of  contact  of  the  cylinder  and  rigid  blo«.k, 
a large  central  portion  of  the  cylinder  »j!J  stari  to  flow 
plastically. 
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